//  only for specific private communication
//
//  Xianguo Lu 
//  lu@physi.uni-heidelberg.de
//  Xianguo.Lu@cern.ch
//

//how to compile
/*
rm -f FitTRDdEdx_cxx.so libFitTRDdEdx.so

aliroot -l -b 
gROOT->LoadMacro("FitTRDdEdx.cxx++")
.q

ln -s FitTRDdEdx_cxx.so libFitTRDdEdx.so
*/

#include "TDirectory.h"
#include "TF1.h"
#include "TH1D.h"
#include "TMath.h"
#include "TMinuit.h"

#include "FitTRDdEdx.h"

#define EPSILON 1e-12

Int_t FitTRDdEdx::fgFitNData=0;
Double_t * FitTRDdEdx::fgFitX = 0x0;
Double_t * FitTRDdEdx::fgFitY = 0x0;
Double_t * FitTRDdEdx::fgFitEy = 0x0;
FitTRDdEdx::FFunc FitTRDdEdx::fgFitFunc = 0x0;

TH1D* FitTRDdEdx::DrawFunc(const TH1 * hh, FFunc func, const Double_t *par, const TString opt, const Int_t color)
{
  TString hname(Form("%sfit",hh->GetName()));
  TH1D *hfit = 0x0;
  gDirectory->GetObject(hname, hfit);
  delete hfit;

  hfit = (TH1D*) hh->Clone( hname );
  for(Int_t ib=1; ib<=hfit->GetNbinsX(); ib++){
    const Double_t ix = hfit->GetBinCenter(ib);
    const Double_t iy = func(&ix, par);
    hfit->SetBinContent(ib, iy);
  }

  hfit->SetLineWidth(2);
  hfit->SetLineColor(color);
  hfit->Draw(opt+" same");

  return hfit;
}

Double_t FitTRDdEdx::ALandau(const Double_t *xx, const Double_t *par)
{
  //
  //p0*landau, npar = 3
  //
  return par[0]*TMath::Landau(xx[0], par[1], par[2], kTRUE);
}

Double_t FitTRDdEdx::ALG(const Double_t * xx, const Double_t * par)
{
  //
  //npar = 4
  //
  const Double_t eps = 1e-6;

  const Int_t ng = 5;
  const Double_t xmin = TMath::Max(0.0, xx[0]-ng*par[3]);
  const Double_t xmax = xx[0]+ng*par[3];

  TF1 ftmp("alg",Form("[0]*TMath::Landau(x, [1], [2], 1)*TMath::Gaus(x, %.10f, [3], 1)", xx[0]), xmin, xmax);

  const Double_t sum = ftmp.Integral(xmin, xmax, par, eps);
  
  return sum;
}

Double_t FitTRDdEdx::NorDimlessLGE(const Double_t * par)
{
  const Double_t xmin = 0;
  const Double_t xmax = 100;

  const Double_t eps = 1e-6;

  TF1 ftmp1("nordimndl1","TMath::Landau(  x, [0], [1]*[0], 1)*TMath::Exp(-[3]*    x/[0])",     xmin,   xmax);
  TF1 ftmp2("nordimndl2","TMath::Landau(1/x, [0], [1]*[0], 1)*TMath::Exp(-[3]*(1/x)/[0])/x/x",    0, 1/xmax);

  const Double_t nor1 = ftmp1.Integral(xmin, xmax, par, eps);
  const Double_t nor2 = ftmp2.Integral(0, 1/xmax, par, eps);
  return nor1+nor2;
}

Double_t FitTRDdEdx::DimlessALGE(const Double_t * xx, const Double_t * par)
{
  const Int_t npar = 5;

  const Double_t eps = 1e-6;
  Double_t absp[npar];
  for(Int_t ii=0; ii<npar; ii++)
    absp[ii] = TMath::Abs(par[ii]);

  const Double_t gs = absp[3]*absp[1];
  const Int_t ng = 5;
  const Double_t xmin = TMath::Max(0.0, xx[0]-ng*gs);
  const Double_t xmax = xx[0]+ng*gs;

  TF1 ftmp("dimlessalgefunc",Form("[0]*TMath::Landau(x, [1], [2]*[1], 1)*TMath::Gaus(x, %.10f, [3]*[1], 1)*TMath::Exp(-[4]*x/[1])", xx[0]), xmin, xmax);

  const Double_t sum = ftmp.Integral(xmin, xmax, absp, eps);
  const Double_t nor = NorDimlessLGE(&(absp[1]));

  return sum/(nor+EPSILON);
}


Int_t FitTRDdEdx::ChisquareFit(const TH1 *hh, FFunc func, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres)
{
  const Int_t nmax = hh->GetNbinsX();
  fgFitX = new Double_t[nmax];
  fgFitY = new Double_t[nmax];
  fgFitEy = new Double_t[nmax];
  fgFitNData = 0;

  for(Int_t ii=hh->GetXaxis()->GetFirst(); ii<=hh->GetXaxis()->GetLast(); ii++){
    const Double_t ey = hh->GetBinError(ii);
    const Double_t yy = hh->GetBinContent(ii);
    //skip empty bin
    if(ey<EPSILON){
      if(yy>EPSILON){
        printf("FitTRDdEdx::ChisquareFit error! %d %e %e\n", ii, ey, yy); exit(1);
      }
      continue;
    }

    //skip low statstics bin
    if(yy<thres){
      continue;
    }

    const Double_t xx = hh->GetBinCenter(ii);

    fgFitX[fgFitNData]=xx;
    fgFitY[fgFitNData]=yy;
    fgFitEy[fgFitNData]=ey;
    fgFitNData++;
  }

  fgFitFunc = func;

  TMinuit * mnt = new TMinuit(npar);
  mnt->SetFCN(FitTRDdEdx::MinFCNChisquare);
    
  const Double_t errdef = 1;
  mnt->SetErrorDef(errdef);
    
  for(Int_t i=0;i<npar;i++){
    const Double_t start=pars[i];
    const Double_t ss=0.01;
    const Double_t low=0;
    const Double_t hi=0;
    
    mnt->DefineParameter(i, Form("p%d",i), start, ss, low,hi);
  }

  const Int_t maxiter = 50000;  
  const Double_t tol= 1;
  Double_t arg[2]={maxiter, errdef*tol};
  Int_t migradflag = -999;
  mnt->mnexcm("MIGRAD", arg, 2, migradflag);

  //get par right after MIGRAD before HESSE modifies them
  Double_t migp[10];
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<npar;ipar++){
    mnt->GetParameter(ipar,migp[ipar], dummy);
  }

  Int_t hesseflag = -999;
  mnt->mnexcm("HESSE", arg ,1, hesseflag);

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("FitTRDdEdx::FitKernel fitting error !! %d %d setting par and err to -999!!!\n", migradflag, hesseflag);
    for(Int_t ipar=0;ipar<npar;ipar++){
      pars[ipar]=-999;
      errs[ipar]=-999;
    }
    chi[0]=chi[1]=-999;
  }
  else{
    for(Int_t ipar=0;ipar<npar;ipar++){
      mnt->GetParameter(ipar, dummy,errs[ipar]);
      pars[ipar]=migp[ipar];
    }
    chi[0] = GetFCNChisquare(npar, pars);
    chi[1] = -mnt->GetNumFreePars();//need to add nx outside FitKernel
  }

  delete mnt;
  delete fgFitX;
  delete fgFitY;
  delete fgFitEy;

  return kfail;
}

Double_t FitTRDdEdx::GetFCNChisquare(const Double_t npar, const Double_t pars[])
{
  Double_t f = 0;

  for(Int_t idata=0; idata<fgFitNData; idata++){
    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t dy = fgFitY[idata] - evaly;
    
    f += dy*dy/(fgFitEy[idata]*fgFitEy[idata]);
  }

  return f;
}

void FitTRDdEdx::MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  f = GetFCNChisquare(npar, pars);
}
